Renta por municipio

0.1 Renta media por municipios

Esta sección presenta un caso de uso en el que aprenderemos a realizar las siguientes tareas básicas:

Importar datos tabulares y datos espaciales.
Realizar un tratamiento de limpieza de datos y cruzar tablas.
Hacer mapas temáticos. Aprenderemos también algunas nociones básicas sobre cómo crear diferentes clases para un conjunto de datos continuo.

Para ello, partiremos de dos ficheros:

Fichero renta_municipio.csv: Este fichero contiene información de la Renta Neta per cápita por municipios, distritos y secciones censales. Esta información se ha extraído del Atlas de distribución de renta de los hogares proporcionado por el INE, y ha sido tratado previamente para adaptar la información al presente ejercicio.
Fichero municipios.gpkg: Es un fichero que contiene datos espaciales (polígonos) de los municipios en España en el año 2019. Se ha extraído del Instituto Geográfico Nacional (IGN) usando el paquete mapSpain.

El primer paso en cualquier tipo de análisis de datos es importar los datos al software de tratamiento (en nuestro caso, R) y analizarlos para conocer el tipo de información que contiene:


# Usaremos paquetes del tidyverse
library(dplyr)
library(readr)


renta <- read_csv("data/renta_municipio.csv", na = ".")

head(renta)
#> # A tibble: 6 × 6
#>   Unidad                          `2019` `2018` `2017` `2016` `2015`
#>   <chr>                            <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
#> 1 44001 Ababuj                        NA     NA     NA     NA     NA
#> 2 4400101 Ababuj distrito 01          NA     NA     NA     NA     NA
#> 3 4400101001 Ababuj sección 01001     NA     NA     NA     NA     NA
#> 4 40001 Abades                     11429  10731  10314   9816   9904
#> 5 4000101 Abades distrito 01       11429  10731  10314   9816   9904
#> 6 4000101001 Abades sección 01001  11429  10731  10314   9816   9904

Podemos comprobar que tenemos información para los años 2015 a 2019. Además, la columna Unidad contiene un literal con el municipio o sección correspondiente.

library(sf)

munis <- st_read("data/municipios.gpkg", quiet = TRUE)

head(munis)
#> Simple feature collection with 6 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -3.140179 ymin: 36.73817 xmax: -2.057058 ymax: 37.54579
#> Geodetic CRS:  ETRS89
#>   codauto ine.ccaa.name cpro ine.prov.name cmun      name LAU_CODE
#> 1      01     Andalucía   04       Almería  001      Abla    04001
#> 2      01     Andalucía   04       Almería  002  Abrucena    04002
#> 3      01     Andalucía   04       Almería  003      Adra    04003
#> 4      01     Andalucía   04       Almería  004 Albanchez    04004
#> 5      01     Andalucía   04       Almería  005 Alboloduy    04005
#> 6      01     Andalucía   04       Almería  006     Albox    04006
#>                             geom
#> 1 MULTIPOLYGON (((-2.775594 3...
#> 2 MULTIPOLYGON (((-2.787566 3...
#> 3 MULTIPOLYGON (((-3.051988 3...
#> 4 MULTIPOLYGON (((-2.181086 3...
#> 5 MULTIPOLYGON (((-2.572442 3...
#> 6 MULTIPOLYGON (((-2.128106 3...

Podemos comprobar que munis es un objeto que contiene Polígonos y varias columnas, entre ellas dos especialmente relevantes: cpro y cmun, que corresponden a los códigos de provincia y de municipio respectivamente. Podemos comprobar que este código también se encuentra en nuestro dataset renta:

# Miro un municipio: Noblejas

renta[grep("Noblejas", renta$Unidad), ]
#> # A tibble: 5 × 6
#>   Unidad                            `2019` `2018` `2017` `2016` `2015`
#>   <chr>                              <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
#> 1 45115 Noblejas                     10591  10314   9751   9484   9124
#> 2 4511501 Noblejas distrito 01       11039  10717  10135   9711   9386
#> 3 4511501001 Noblejas sección 01001  11039  10717  10135   9711   9386
#> 4 4511502 Noblejas distrito 02       10276  10029   9475   9319   8938
#> 5 4511502001 Noblejas sección 02001  10276  10029   9475   9319   8938

munis[grep("Noblejas", munis$name), c("name", "cpro", "cmun")]
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -3.489824 ymin: 39.93003 xmax: -3.372611 ymax: 40.05017
#> Geodetic CRS:  ETRS89
#>          name cpro cmun                           geom
#> 4985 Noblejas   45  115 MULTIPOLYGON (((-3.44681 40...

En el caso de Noblejas, el código completo es 45115. Sin embargo, en el caso de la tabla renta, debemos extraer ese valor del literal. Para ello debemos manipular la columna y extraer la primera palabra de la columna Unidad:


# Creo una función y la aplico a toda la columna
extrae_codigo <- function(x) {
  unlist(strsplit(x, " "))[1]
}

renta$codigo_ine <- sapply(as.character(renta$Unidad), extrae_codigo)

head(renta[c("Unidad", "codigo_ine")])
#> # A tibble: 6 × 2
#>   Unidad                          codigo_ine
#>   <chr>                           <chr>     
#> 1 44001 Ababuj                    44001     
#> 2 4400101 Ababuj distrito 01      4400101   
#> 3 4400101001 Ababuj sección 01001 4400101001
#> 4 40001 Abades                    40001     
#> 5 4000101 Abades distrito 01      4000101   
#> 6 4000101001 Abades sección 01001 4000101001

Ahora, es necesario crear la misma variable en munis para poder realizar el cruce:


munis$codigo_ine <- paste0(munis$cpro, munis$cmun)

head(munis[, c("name", "codigo_ine")])
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -3.140179 ymin: 36.73817 xmax: -2.057058 ymax: 37.54579
#> Geodetic CRS:  ETRS89
#>        name codigo_ine                           geom
#> 1      Abla      04001 MULTIPOLYGON (((-2.775594 3...
#> 2  Abrucena      04002 MULTIPOLYGON (((-2.787566 3...
#> 3      Adra      04003 MULTIPOLYGON (((-3.051988 3...
#> 4 Albanchez      04004 MULTIPOLYGON (((-2.181086 3...
#> 5 Alboloduy      04005 MULTIPOLYGON (((-2.572442 3...
#> 6     Albox      04006 MULTIPOLYGON (((-2.128106 3...

Ya estamos listos para realizar el cruce. Además, seleccionaremos sólo las columnas que vamos a usar, en este caso la del año 2019:


munis_renta <- munis %>%
  left_join(renta) %>%
  select(name, cpro, cmun, `2019`)

Cuando crucemos datos espaciales con datos no espaciales en R, es importante que el primer dataset sea el que contiene los datos espaciales. Esto es así porque el objeto resultante “hereda” la clase del primer objeto. A modo de ejemplo, si realizáramos el proceso poniendo los datos espaciales en el lado derecho del join, los datos finales no serán espaciales:


# Miramos la clase de munis_renta

class(munis_renta)
#> [1] "sf"         "data.frame"

# Es un sf, por tanto espacial

# ¿Que pasa si realizamos el cruce de la otra manera?
renta %>%
  left_join(munis) %>%
  select(name, cpro, cmun, `2019`) %>%
  class()
#> [1] "tbl_df"     "tbl"        "data.frame"

# Es un tibble o data.frame, pero no es espacial!

Una vez que tenemos los datos unidos podemos realizar algunos análisis básicos, como la realización de un histograma


library(ggplot2)

munis_renta %>%
  ggplot(aes(x = `2019`)) +
  geom_histogram(color = "darkblue", fill = "lightblue") +
  scale_x_continuous(labels = scales::label_number_auto()) +
  scale_y_continuous(labels = scales::label_percent()) +
  labs(
    y = "",
    x = "Renta neta media por persona (€)"
  )

Podemos observar que la renta presenta una distribución Gamma con un gran de municipios concentrados en zonas medias de renta y pocos municipios en tramos de rentas altas. Como veremos más adelante, esta distribución va a afectar a la información que transmite el mapa.

Vamos a realizar ahora un mapa de coropletas mostrando la distribución de la renta usando los valores brutos de renta sin modificar:


ggplot(munis_renta) +
  # Usamos geom_sf, y como aes() lo que queremos mostrar, en este caso, el
  # color del polígono representa la renta. Vamos a retirar los bordes con
  # color = NA
  geom_sf(aes(fill = `2019`), color = NA) +
  theme_minimal() +
  scale_fill_continuous(labels = scales::label_number(
    big.mark = ".",
    decimal.mark = ",",
    suffix = " €"
  )) +
  labs(
    title = "Renta neta media por persona",
    caption = "Datos: INE"
  )

Este primer mapa no es demasiado informativo, por los siguientes motivos:

Existe una serie de municipios para los que no tenemos datos.
La escala de color no es la más adecuada.
Dada la distribución de los datos, puede ser adecuado crear grupos de renta para que el mapa sea más interpretable.

Vamos a probar a eliminar los municipios sin datos y a cambiar la escala de color:


munis_renta_clean <- munis_renta %>% filter(!is.na(`2019`))

ggplot(munis_renta_clean) +
  geom_sf(aes(fill = `2019`), color = NA) +
  # Cambiamos la paleta de colores, vamos a usar una paleta denominada Inferno,
  # ya incluida en base R con hcl.colors

  # Como son datos continuos, puedo usar Inferno
  scale_fill_gradientn(
    colours = hcl.colors(20, "Inferno", rev = TRUE),
    labels = scales::label_number(
      big.mark = ".",
      decimal.mark = ",",
      suffix = " €"
    )
  ) +
  theme_minimal() +
  labs(
    title = "Renta neta media por persona",
    caption = "Datos: INE"
  )

Este mapa nos da algo más de información, y parece intuirse que las rentas más altas se encuentran en zonas de País Vasco, Madrid y Cataluña. Sin embargo, el hecho de que la distribución de los datos no sea normal está afectando a la visualización.

Para intentar atajar este problema, podemos dividir nuestros datos en clases, por ejemplo cuartiles o deciles. Existen varios métodos de clasificación de datos, que en R se encuentran implementados en el paquete classInt. A continuación vamos a plantear diversos métodos de clasificación y observar cómo la “historia” que cuenta el mapa varía en función de dichas clases. En este ejemplo planteamos los siguientes métodos de clasificación:

El método de deciles: Consiste en crear 10 categorías incluyendo el mismo número de registros en cada una de ellas.
El método de intervalos equivalentes: Este método divide el rango de valores en un número de grupos definido. La distancia de todos los intervalos es idéntica, por lo que este método no tiene en cuenta la distribución de los registros.
El método de Fisher-Jenks: Este método se desarrolló específicamente para la clasificación de datos espaciales y su visualización en mapas. Produce agrupaciones de tal manera que los datos de cada grupo son “cercanas” entre sí y sustancialmente distintas de los valores de otros grupos.


library(classInt)
# Vamos a probar 3 métodos de clasificación: Deciles, tramos de Renta
# equidistantes y Fisher and Jenks

deciles <- classIntervals(munis_renta_clean$`2019`,
  style = "quantile", n = 10
)
deciles
#> style: quantile
#>     [5898,8935.6)   [8935.6,9662.2)  [9662.2,10352.8)   [10352.8,10918) 
#>               656               656               655               654 
#>     [10918,11462)   [11462,11998.6) [11998.6,12651.4) [12651.4,13475.8) 
#>               655               658               656               655 
#> [13475.8,14618.4)   [14618.4,26367] 
#>               656               656
plot(deciles, pal = hcl.colors(20, "Inferno"), main = "Deciles")



# Tramos equidistantes en términos de renta
equal <- classIntervals(munis_renta_clean$`2019`,
  style = "equal", n = 10
)
equal
#> style: equal
#>     [5898,7944.9)   [7944.9,9991.8)  [9991.8,12038.7) [12038.7,14085.6) 
#>               103              1510              2374              1637 
#> [14085.6,16132.5) [16132.5,18179.4) [18179.4,20226.3) [20226.3,22273.2) 
#>               702               161                52                14 
#> [22273.2,24320.1)   [24320.1,26367] 
#>                 3                 1
plot(equal, pal = hcl.colors(20, "Inferno"), main = "Equidistantes")


fisher <- classIntervals(munis_renta_clean$`2019`,
  style = "fisher",
  # Fuerzo para mejorar la comparación entre métodos
  n = 10
)
fisher
#> style: fisher
#>     [5898,8782.5)   [8782.5,9806.5)    [9806.5,10755)     [10755,11692) 
#>               530               911               974              1158 
#>   [11692,12738.5) [12738.5,13909.5)   [13909.5,15257)   [15257,17188.5) 
#>              1101               868               595               298 
#> [17188.5,20894.5)   [20894.5,26367] 
#>               111                11
plot(fisher,
  pal = hcl.colors(20, "Inferno"),
  main = "Fisher-Jenks"
)

Podemos observar lo siguiente:

El último decil de renta se corresponde a un rango de entre 15.000 y 25.000 €.
El método por deciles proporciona unos grupos con valores de renta muy parecidos entre sí en los valores medios. Esto es debido a la propia distribución de la variable.
El método de rangos equidistantes proporciona algunos grupos con un número muy reducido de municipios.
El método de Fisher-Jenks puede proporcionar unas clases con unos rangos más apropiados para los tramos altos de renta.

Vamos ahora a realizar 3 mapas distintos, creando clases de renta según cada uno de los métodos anteriormente mostrados.

Deciles

# Extraigo los valores de corte
breaks_d <- deciles$brks

# Y creo unas etiquetas básicas para cada clase
# Creo una función específica para crear etiquetas formateadas
label_fun <- function(x) {
  l <- length(x)
  eur <- paste0(prettyNum(round(x, 0),
    decimal.mark = ",",
    big.mark = "."
  ), " €")

  labels <- paste(eur[-l], "-", eur[-1])
  labels[1] <- paste("<", eur[1])
  labels[l - 1] <- paste(">", eur[l - 1])
  return(labels)
}

labels_d <- label_fun(breaks_d)

munis_renta_clean$Deciles <- cut(munis_renta_clean$`2019`,
  breaks = breaks_d,
  labels = labels_d,
  include.lowest = TRUE
)

ggplot(munis_renta_clean) +
  # Cambiamos la variable que usamos para crear el mapa
  geom_sf(aes(fill = Deciles), color = NA) +
  # Necesito cambiar el scale, ya no es continua
  scale_fill_manual(values = hcl.colors(length(labels_d),
    "Inferno",
    rev = TRUE
  )) +
  theme_minimal() +
  labs(
    title = "Renta neta media por persona",
    caption = "Datos: INE"
  )

Este mapa ya nos permite observar patrones geográficos, donde se ve una clara diferencia entre la Comunidades Autónomas del Norte y las del Sur. Veamos una representación distina usando otras clases diferentes:


breaks_e <- equal$brks
labels_e <- label_fun(breaks_e)

munis_renta_clean$Equal <- cut(munis_renta_clean$`2019`,
  breaks = breaks_e,
  labels = labels_e,
  include.lowest = TRUE
)

ggplot(munis_renta_clean) +
  # Cambiamos la variable que usamos para crear el mapa
  geom_sf(aes(fill = Equal), color = NA) +
  scale_fill_manual(values = hcl.colors(length(labels_e),
    "Inferno",
    rev = TRUE
  )) +
  theme_minimal() +
  labs(
    title = "Renta neta media por persona",
    caption = "Datos: INE"
  )

Este otro mapa, sin embargo, se parece más al mapa que hicimos con los datos sin clasificar, donde el peso visual se concentra más bien en los municipios con rentas mucho más altas que el resto (por encimoa de 18.000 €).

Veamos el mismo mapa usando Fisher-Jenks:


breaks_f <- fisher$brks
labels_f <- label_fun(breaks_f)

munis_renta_clean$`Fisher-Jenks` <- cut(munis_renta_clean$`2019`,
  breaks = breaks_f,
  labels = labels_f,
  include.lowest = TRUE
)

ggplot(munis_renta_clean) +
  # Cambiamos la variable que usamos para crear el mapa
  geom_sf(aes(fill = `Fisher-Jenks`), color = NA) +
  scale_fill_manual(values = hcl.colors(length(labels_f),
    "Inferno",
    rev = TRUE
  )) +
  theme_minimal() +
  labs(
    title = "Renta neta media por persona",
    caption = "Datos: INE"
  )

En este mapa se puede observar de una manera más clara un cluster adicional de renta en la zona de Asturias y el norte de León. Además, gracias a la escala de colores puede intuirse que este clúster de renta no presenta valores tan altos como los observados en País Vasco o Madrid.

En conclusión, en el momento de realizar una visualización de datos es importante conocer el dato a representar, así como entender algunas propiedades básicas de la distribución subyacente. También hemos podido observar que hay ciertas decisiones estéticas (datos continuos vs. agrupados, escala de colores) que tienen una influencia significativa en cómo se percibe la información representada. Es responsabilidad del creador de la visualización el conocer todos estos factores y aplicarlos de manera conveniente.

Como ejemplo, aquí tenemos un mapa de los mismos datos pero la información se presenta de manera sesgada, ¿puedes identificar los motivos?

Ejemplo de visualización sesgada

En este mapa parece que la renta per cápita de las comunidades del norte es desproporcionadamente superior a las del sur. Resumimos aquí los sesgos introducidos en el mapa:

En primer lugar, se han creado un elevado número de grupos en las zonas de rentas bajas. De esta manera, la escala del mapa parece estar muy fragmentada, sin embargo muchos de esos grupos apenas contienen municipios. A modo de ejemplo, los primeros cuatro grupos únicamente contienen 32 municipios.
Los grupos no se adaptan a la distribución subyacente de los datos. La mediana de los datos (11.462 €) estaría situada en la antepenúltima clase, de manera que los dos grupos de mayor renta contienen el 50% de los municipios.
Además, la escala de color se ha manipulado, de manera que los grupos de mayores renta destaquen más que el resto de manera notoria.

Renta por municipio

Diego Hernangómez

2022-01-14

0.1 Renta media por municipios